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Grains inside a vertically vibrated box undergo a transition from a density inverted and hor¬ 
izontally homogeneous state, referred to as the granular Leidenfrost state, to a buoyancy-driven 
convective state. We perform a simulational study of the precursors of such a transition, and 
quantify their dynamics as the bed of grains is progressively fluidized. The transition is preceded 
by transient convective states, which increase their correlation time as the transition point is ap¬ 
proached. Increasingly correlated convective flows lead to density fluctuations, as quantified by the 
structure factor, that also shows critical behaviour near the transition point. The amplitude of the 
modulations in the vertical velocity field are seen to be best described by a quintic supercritical 
amplitude equation with an additive noise term. The validity of such an amplitude equation, and 
previously observed collective semi-periodic oscillations of the bed of grains, suggests a new inter¬ 
pretation of the transition analogous to a coupled chain of vertically vibrated damped oscillators. 
Increasing the size of the container shows metastability of convective states, as well as an overall 
invariant critical behaviour close to the transition. 

PACS numbers: 45.70.Qj, 46.65.+g, 05.40.-a 


I. INTRODUCTION 


Granular materials, defined as collections of dissipa¬ 
tive particles large enough so that thermal fluctuations 
can be ignored, are an archetypal study case of complex 
dynamical systems. Decades of research have revealed 
many novel non-equilibrium phase transitions and col¬ 
lective behaviors m, the study of which not only has 
a fundamental physical interest, but is also relevant for 
many industries [THg. Many of these behaviors show a 
striking similarity with molecular fluids or solid phenom¬ 
ena HDHIS], and some have even been successfully de¬ 
scribed by equilibrium theories M- Studying the origin 
of these agreements advances our understanding of far- 
from-equilibrium states, and explores the limits of contin¬ 
uum descriptions of discrete systems. Furthermore, the 
low number of constituents, when compared to molecu¬ 
lar counterparts, makes granular materials particularly 
suited for the study of noise effects in spatially extended 
transitions, a subject of increasing physical interest due 
to the ubiquitous presence of fluctuations in natural phe¬ 
nomena [smSHIH]. 

In order to keep granular media fluidized it is necessary 
to provide energy to the system. Previously this has been 
done in several distinct ways, as for example electromag- 
netically HHHD], by shearing |21j , or by boundary forces 
such as rotating a drum |22j or vibrating the grains’ con¬ 
tainer |2fl] . In vertically vibrated systems several com¬ 
plex collective dynamic behaviors have been observed, 
such as segregation [53J [^, pattern formation [B] and 
phase separation [5]. One particular case of the latter is 
the granular Leidenfrost state, where a dense, solid- or 
fluid-like region is sustained by a highly agitated low den¬ 
sity gaseous region in contact with the vibrated bottom 


wall [211123 ■ It is so called because of the clear anal¬ 
ogy with the water-over-vapor phenomenon observed in 
molecular fluids in contact with a high temperature sur¬ 
face [28]. If the vibration strength is increased, the Lei¬ 
denfrost state evolves to a buoyancy-driven convective 
state [21]) in analogy to Rayleigh-Bernard convection. 

In the following work we study the precursors of the 
transition from the granular Leidenfrost to the buoyancy- 
driven convective state in the context of bifurcations and 
critical theory. Previous experimental and simulational 
works determined the transition points as a function of 
the energy injection and the amount of particles in the 
system [23 ES]- H also shown that granular hydrody¬ 
namics is able to quantitatively capture the critical points 
of this instability, by performing a linear stability analysis 
of perturbations over the Leidenfrost state [221 [30] . Here 
we explore further the regions close to the transition, 
motivated by the presence of complex transient dynamics 
which are expected to be dominated by fluctuations. This 
transition is an excellent candidate for studying the in¬ 
fluence of fluctuations in hydrodynamic-like instabilities, 
due in part to its similarity with the Rayleigh-Benard 
instability present in regular fluids. Our goal is to in¬ 
crease the knowledge about the origin and evolution of 
the perturbations that lead to the instability, from both 
the microscopic and macroscopic perspectives, and re¬ 
late the transition to other analogous dynamics through 
an unstable-mode amplitude equation. 


After specifying the system and simulation methods 
(Sec. we begin by characte rizing the two states in¬ 
volved in the transition (Sec. Ill A I, and determining 
the phase space of the system by means of a convec¬ 
tion intensity order parameter (Sec. IIIB). With this, we 
are then able to study time-dependent transient convec- 
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tive states, that are present far below the transition and 
show a critically increasing correlation time as the tran¬ 
sition is approached (Sec. IIIC). Furthermore, the static 
structure function allows us to study the evolution of the 
relevant length-scale in this pattern formation scenario, 
and see its behaviour prior to the transition (Sec. |IIID ). 
Finally, it is observed that the amplitude of the criti¬ 
cal pattern follows a growth ratio that is consistent with 
a quintic supercritical bifurcation, associated with para¬ 
metrically driven spatially extended systems [niEiiEg 
(Sec. IIIF). We suggest that the agreement with this 
universality class comes from the presence of collective 
semi-periodic oscillations, so called low-frequency oscil¬ 
lations (LFOs), present in density inverted systems [33]. 
All results are presented for different boundary condi¬ 
tions and sizes of the container, allowing us to observe 
the influence of confinement and variations of the total 
number of particles. 


narrow 



FIG. 1: Schematic representation (not to scale) of the setup. 
Two different geometries are considered: narrow (top) and 
wide (bottom). Lengths are given in units of particle diame¬ 
ters d. 


II. SYSTEM AND SIMULATIONS 

The setup consists of a quasi-two-dimensional rectan¬ 
gular box with open top, vibrated in the vertical direc¬ 
tion. Two different box widths are considered, defining 
the narrow system, with = 50, and the wide one, with 
lx = 400. The depth of the container, on the other hand, 
is kept constant, ly = 5] a, schematic representation of the 
studied geometries is shown in Fig. Here, and in what 
follows, we use dimensionless quantities with d as length- 
scale and y/djg as timescale, and thus \fgjd as velocity 
units; when necessary, dimensional quantities will be dis¬ 
tinguished by a tilde, i.e. lx = Ixd. Grains are considered 
to be perfectly spherical, frictionless and monodisperse in 
size and mass. Their total number N is determined by 
the number of filling layers F = N/{lxly), which we fix at 
F = 12. Previous studies show that both the Leidenfrost 
and the buoyancy-driven convective states are observable 
for this number of layers [M] . The whole box (base and 
side walls) is vertically vibrated in a bi-parabolic, quasi- 
sinusoidal way with a given frequency uj and amplitude 
A. The use of a quadratic interpolation instead of a sine 
function gives a considerable speed advantage in simula¬ 
tions, as the collision times with the moving walls can be 
predicted analytically. Previously, test simulations have 
been done using a sine function for exemplary cases, and 
no significant difference was observed [33 ES] . The am¬ 
plitude of oscillation is kept fixed, A = 0.1, and thus the 
energy injection is controlled by the angular frequency uj. 
The low amplitude is chosen to reduce as much as possi¬ 
ble geometrical effects of the moving boundary (such as 
shock waves [36] L and approximate the limit of a temper¬ 
ature boundary condition m- Moreover, low amplitudes 
eliminate other inhomogeneous states for lower energies, 
such as undulations [34], which are not the object of this 
study. Overall, our selection of parameters is based on 
previous experimental setups where the transition was 
previously reported [13 El] ■ 


The system is simulated using an event-driven (ED) 
hard-sphere algorithm. The advantage of using ED sim¬ 
ulations over regular time-stepping methods is straight¬ 
forward: computational speed. Even though the number 
of particles is relatively low (^ 10^), the high frequen¬ 
cies and very long physical times, of the order of hours, 
make the use of discrete particle methods (DPM) infea¬ 
sible. In DPM simulations time-steps are constant and 
should be at least one order of magnitude lower than the 
collision time, which in itself must be at least an order of 
magnitude smaller than the lowest relevant time-scale, in 
our case T = 2tt/uj [38]. Thus, for the high frequencies 
considered in our study, the small time-step prohibits to 
simulate in a practical time the long transients involved 
near a transition. On the contrary, the average time-step 
in ED is determined mainly by the density of the system, 
and not directly dependent on the frequency of oscillation 
of the container. 

Collisions between particles are modeled by a normal 
restitution coefficient, rp = 0.9 [39]. In order to avoid 
inelastic collapse, the TC model is used, where particle 
collisions are considered elastic if they occur within a 
given time, which we take as tc = 10“® [40]. This essen¬ 
tially sets a lower limit for physically relevant velocities, 
as also slightly decreases the packing fraction of high den¬ 
sity regions; possible relevant effects will be noted when 
appropriate. 

Regarding boundary conditions, we consider both 
cases of periodic (PBC) and solid boundary conditions, 
with either elastic or dissipative walls (EEC and DBG, 
respectively). The different boundary types are only 
applied in the ^-direction, as we would like to investi¬ 
gate the effects they have on the transition independent 
of other factors, as increased overall dissipation or free- 
volume; setting dissipative or periodic boundaries also in 
the y-direction would make the comparison less straight¬ 
forward. Dissipative walls are set with the same resti- 
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FIG. 2: In the narrow system, time averaged number density 
of particles {n{x,z))^ (fop)? granular temperature {T{x,z))^ 
(middle) and velocity field {v{x,z))^ (bottom), for systems 
in the granular Leidenfrost state (left) and in the buoyancy- 
driven convective state (right). 


ing velocity, SksT = The fields clearly 

show that in the Leidenfrost state, both p and T are 
homogeneous in the :c-direction, while in the convective 
state the profiles are modulated by a dominant mode kc- 
That is, the transition is morphogenetic |42] . as a pattern 
or new relevant length-scale arises from a homogeneous 
state. The convective mode defines the typical size of a 
pair of convective cells, Ac; in the case shown in Fig. 
{Xc)t ~ 50, that is, kc = 1/Ac ~ 0.02. 

It is important to remark that the buoyancy-driven 
convective state is also density inverted (see Fig. If), 
and thus this characteristic is not a sufficient condition 
to define the Leidenfrost state. We demand two further 
properties for the system to be considered in this state: 
(a) higher density regions present distinct dynamics to 
the lower density ones (gas/fluid or gas/solid), to distin¬ 
guish it from completely gaseous states [43]; and (b) the 
system remains horizontally homogeneous, to differenti¬ 
ate it from the convective state. In short, we define the 
granular Leidenfrost state as a density inverted, phase 
coexisting, horizontally homogeneous state. 

As the energy input increases, the bed of grains in 
the dense region progressively looses its horizontal homo¬ 
geneity, giving rise to convection; this is what we refer 
to as the granular Leidenfrost to buoyancy-driven con¬ 
vection transition or, in short, the LBC transition. In 
the following we define an order parameter based on the 
evolution of the velocity field, and observe its behaviour 
through the transition. 


tution coefficient as between particles, = 0.9. The 
effects of dissipative walls on convective states have al¬ 
ready been studied in similar setups, both experimentally 
and numerically |41] . Here we are interested in the effects 
of walls on the excitation or suppression of the modes rel¬ 
evant in the transition. Elastic walls (EEC) are used in 
order to see the influence of excluded volume effects near 
the sidewalls when comparing with periodic walls (PBC), 
as also to facilitate the analysis of fluctuations by fixing 
a reference frame. Furthermore, PBCs are used to study 
the dynamics of the bed of grains without confinement. 


III. RESULTS 
A. Macroscopic description 

The most evident difference between the granular Lei¬ 
denfrost and buoyancy-driven convective states is the 
level of horizontal homogeneity. Fig. shows time av¬ 
eraged number density (n(x,z))j, granular temperature 
(T(x,z))j and velocity helds {v{x,z))^ in each state, for 
narrow systems with EBC. The fields are obtained by 
binning the system in squares of size d. Time averages, 
()(, are always taken for at least 10^^/d/g, which in di¬ 
mensional terms for d = Imm would correspond to ex¬ 
periments of about Hfteen minutes. The granular tem¬ 
perature is defined as the kinetic energy of the fluctuat- 


B. Convection intensity 

For the study of critical behaviors it is of fundamental 
importance that the transition region between the two 
states is accurately measured. The different states can 
be easily distinguished by looking at the time-average 
velocity fields, which suggest the use of the convection 
intensity order parameter, defined as 

C = imax^(max,j(w^(x, z)) - mina;(wz(x, z))). (1) 

Here Vz{x,z) is the scalar field of velocities in the z- 
direction, and the maxima are taken first over z and then 
over X. In words, C corresponds to half the highest dif¬ 
ference of the vertical velocities at a particular height 
of the container. In a convective state C is expected to 
be significantly higher than in a random flux case, due 
to the presence of stable upwards and downwards flux 
regions (as can be seen in Fig. [^). Even though the av¬ 
erage vertical velocity is expected to scale with Auj, the 
localization of the energy fluxes in the convective states 
is what produces a higher deviation, and thus a higher 
C. The time averaged convection intensity, (C')j, cap¬ 
tures the transition as a rapid increase with w, as shown 
in Fig. for all considered systems. 

In the Leidenfrost state (C')^ increases linearly with oj, 
and is lower than the characteristic velocity of energy in¬ 
jection, Auj. This is followed by a transition region, were 
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FIG. 3: Time averaged convection intensity (C')^, defined in 
the main text as a function of the angular driving fre¬ 
quency ui, for the narrow (top) and wide (bottom) containers 
with the boundary conditions indicated in the labels. Dashed 
lines indicate the transition region for the corresponding BC, 
as specified in the main text. The thick gray line corresponds 
to Auj, the characteristic shaking velocity. The insets show 
the convection intensity normalized by the driving frequency, 
C* = (C)^ /Alo, as a function of the bifurcation parameter 
e = {lj — uJc)luJc- 


(C)j increases sharply and superlinear on Auj, eventu¬ 
ally surpassing the Auj line. Finally, (C)^ saturates as 
the system enters the stable buoyancy-driven convective 
state. Quantitatively, we define the limits of the tran¬ 
sition region by looking at the intersection of the ini¬ 
tial and final linear behaviors with the increasing tran¬ 
sient behaviour, the two points defining the width of the 
transition, 5uj, and their average the critical frequency 
Wc, which coincides within measurement error with the 
condition (C')j = Auj. For the narrow container, this 
results in critical frequencies and widths of transition 
Wc = 33.4 ± 0.1, 5uj = 0.22 ± 0.01 for EEC and PEC, 
and Wc = 32.9 ± 0.1, 5uj = 0.26 ± 0.01 for DEC, with 
the error given by the resolution of the simulations in w. 
That is, elastic boundary conditions have no measurable 
infiuence when compared to periodic boundaries, which 
suggests that excluded volume effects due to the presence 
of walls can be disregarded already for Ij, = SOd. Dissi¬ 
pative boundaries, on the other hand, have the at first 
counterintuitive effect of decreasing Wc, while increasing 
Suj; even though overall the system presents more dis¬ 
sipation compared to the EEC case, inelastic sidewalls 
slightly reduce the energy needed to trigger the transi¬ 
tion compared to elastic boundaries. 

Eoundary conditions in the wide container become ir¬ 
relevant, with all cases given by ujc = 33.0 ±0.1 and 
Suj = 0.29 ± 0.02. Quantitatively, the critical points are 
slightly lower and the transitions wider, which we believe 


is due to the infiuence of the confinement in the narrow 
container. It is worthy to remark that the amount of 
energy needed for the creation of the convective cells is 
practically invariant on or, equivalently, the number of 
convective rolls ric = suggesting that the inter¬ 

action between rolls has no influence on their creation. 
Nevertheless, we notice that when EEC or DEC are used, 
convection cells are seen to appear first at the boundaries, 
and the boundary rolls are more stable when compared 
to the bulk of the system. This, nevertheless, happens at 
the same uj^ as with PEC, suggesting that solid bound¬ 
aries have no relevant influence on the flux (nv) strength, 
but do promote the appearance of convective cells near 
them. 


When normalized by Auj, we can recognize in C* = 
(C)j /Auj a shape characteristic of a supercritical pitch- 
fork bifurcation, as shown in the insets of Pig. The 
second branch of the ideal pitchfork supercritical bifur¬ 
cation would correspond to taking the minimum in x, 
instead of the maximum, in Q. When the bifurcation 
parameter e = (w — uJc)/uJc is used as control parame¬ 
ter, all three boundary condition cases coincide for all 
system sizes considered. This suggests that the tran¬ 
sition presents universal behaviour, independent of the 
amount of dissipation. It is also a confirmation that the 
critical points are well defined. With the phase-space 
determined, next we characterize the precursors of the 
transition by looking first at correlations of the velocity 
field (IIIC), and then at den sity fl uctuations by means 
of the static structure factor (HID). 


C. Time-dependent fluctuating convective flows 

Far below the transition point in the Leidenfrost state, 
starting from e > —0.5, time-dependent fluctuating con¬ 
vective flows can be observed. These are analogous to the 
precursor’s fluxes present in the classical fluid Rayleigh- 
Eenard convection transition [44] . which were theoreti¬ 
cally predicted and relatively recently observed by care¬ 
ful experiments in gaseous media |45| . In our case, the 
convective rolls can be easily identified when observ¬ 
ing the evolution of the short-time-averaged transient 
velocity fields, as shown in Fig. The cells are con¬ 
stantly generated anywhere in the container, but more 
frequently next to walls, this of course when they are 
present, i.e. in the EEC and DEC cases. Two fundamen¬ 
tal aspects differentiate such a transient state from the 
fully developed buoyancy-driven convective state above 
e = 0. Pirst, the circulation of particles is not associated 
with mean density or temperature inhomogeneities (it is 
time-dependent). Second, the convective velocity field is 
present only as an average, and thus is not correlated 
with the instantaneous velocity of the particles. That 
is, the velocities of the fluctuating convective flows are 
much smaller than the amplitude of the fluctuating ve¬ 
locities {Vt), in contrast to the buoyancy-driven convec¬ 
tion case, where they are comparable (see Fig. [^. This 
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FIG. 5: (a) Velocity correlation functions for several uj and 
EEC in the narrow container, (b) Characteristic time-scale of 
fluctuating convection Tv, corresponding to the exponent of 
the long term exponential decay of the self correlation func¬ 
tion F^, as a function of the bifurcation parameter e. The 
dashed lines indicate best fits of the form indicated in the 
main text. 


FIG. 4: Transient velocity fields for e = —0.15, each averaged 
over 5 oscillation periods, showing the emergence and decay 
of a fluctuating convective cell in a section of a wide container. 
From blue to red, the color and size of the vectors corresponds 
to their norm. 

has the consequence that, as there is no localization of 
the fluxes, their effect is not reflected in (C')^. 

In order to characterize the stability of the transient 
convective cells, the self-correlation of the fluctuating ve¬ 
locity field is computed, 

Fv{t) = cf {Sv{x, t + r) ■ dv{x, t))^ 

with 5v = v{x,t) — {v{x,t))^^ and cp a normalization 
constant such that Fy{0) = 1. Fig. shows Flu(r) for 
characteristic cases of uj. In the following we focus only 
on EEC and DEC, as they considerably simplify the com¬ 
putation of self-correlation functions by impeding the 
convective rolls to drift in the cc-direction, as they do 
with PEC. Visual inspection and preliminary analysis of 
the PEC case suggest that the results can be generalized 
to this case as well. All correlations present a common 
shape: an initial quick, power-law-like decay followed by 
a slower exponential decrease. The rapid decorrelation at 
short time-scales confirms that the particles’ instant ve¬ 
locities are mostly fluctuating, and do not present a high 
time correlation. On the other hand, for longer times the 
correlation is comparatively lower, but still considerable, 
and decays slower. This is a signal of long-term aver¬ 
age preferred fluxes. As expected by the critical slowing 
down of fluctuations near the transition, the overall cor¬ 
relation of this region increases as the critical point is 
approached, as can also be seen in Fig. [Da- The charac¬ 
teristic time of decorrelation Ty is obtained by consider¬ 
ing Fy - exp(-T/T„). Fig. shows ujTy as a function 
of e, from where we find a powerlaw Ty ^ e~^/uj with 


exponent ^ ^ 0.59 ± 0.02. Closer to the critical point 
the measurement error becomes significant. The data is 
presented for the whole range in uj where the Leidenfrost 
state is present, which is one and a half decades in e. 

Wide systems present the same overall features as the 
narrow container; Ty can be determined with a higher 
precision -as noise is reduced with a higher number of 
particles-, and presents the same, within error, critical 
exponent ^ as in the narrow case, ^ ~ 0.60 ± 0.02. That 
is, transient convective flows are independent of the size 
of the container. 

Also visible in Fi,(r) are wide peaks at regular in¬ 
tervals, signals of a quasi-periodic time-scale of corre¬ 
lation. Ey observing the evolution of the center of mass, 
and computing its fast-Fourier transform, it was verified 
that this periodic correlation corresponds to the recently 
reported low-frequency oscillations, present in density 
inverted agitated systems [33l [46]. The quasi-periodic 
movement is coupled with a breathing behaviour of the 
dense bed of grains, which increases and decreases its 
granular temperature. Here we do not analyze this fur¬ 
ther; for a detailed study of the phenomena we refer the 
reader to |33j . and a further experimental study in |46j . 


D. Static structure function 

As energy input increases, for e > — O.I, density fluc¬ 
tuations arise, clearly recognizable as modulations in the 
surface of the bed of grains. To analyze their behaviour 
we compute the static structure function, 

S{k) = ^ 


( 2 ) 
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FIG. 6: Structure factor, S{k), for narrow (left) and wide 
(right) containers for the bifurcation parameters specified. 
Dashed lines correspond to PBC, while solid lines have EBC. 
The vertical solid line indicates the l/Za, point. 

with h the Fourier components of the depth-averaged 
number density field in the x-direction, 

G/<5x 

n{k,t) = (3) 

3 

Notice that we define k = 1/A, for a more straightforward 
comparison between wave number k and wavelength A. 
The position Xj is given by regular intervals, Xj = ^dx + 
jSx, with Sx = 0.1 the coarse graining length. Notice 
that instead of considering the particles’ position in the 
definition of n we use the averaged density profiles, as 
it significantly increases the speed of computation. This 
approximation holds only for low wave-numbers, that is, 
1/fc ^ 6x, which is the region we are interested in. Test 
cases were done with the usual definition with particle 
positions, and no significant differences were observed. 

Transient modulations of the bed are captured in S{k) 
by the appearance and steady increase of a narrow peak 
at /c « 0.02, as shown in Fig. |^for both the narrow and 
wide containers. We define the critical mode kc by the 
position of this maximum, that is Sm = max(5'(fc)) = 
S'(fcc). Thus, the associated wavelength at the transition 
point. A* « 50, corresponds to the size of the smallest 
stable convection roll, seen to be independent of lx for 
lx > X*. Notice that this corresponds to Uc = 2 for the 
narrow container, and ric = 8 for the wide container. 
What other factors may affect A* is not studied further 
here, although we notice that previously realized stabil¬ 
ity analysis of the granular hydrodynamic equations have 
found an expression for A* as a function of the constitu¬ 
tive relations, which are in themselves dependent on the 
particle properties [3Dj. 

Notice from Fig. |^that for e = —0.1 the correlation of 
the transient convective flows was significant, but S{k) 
has no relevant maximum. This confirms that fluctuat¬ 
ing convective flows take place in a stable homogeneous 
Leidenfrost state, and are not accompanied by any rele¬ 
vant excitation of the critical mode in the density (and 
temperature) field. 


FIG. 7: The most unstable mode kc, defined by the maximum 
of the structure factor max(S(k)) = S{kc), as a function of 
the bifurcation parameter, for narrow (left) and wide (right) 
containers and the boundary conditions specified. The in¬ 
set shows the convective length-scale Ac as a function of the 
shaking strength, as defined in the main text. 

Previous simulational and experimental works have 
stated that Ac scales linearly with the shaking strength 
S = A^Co'^/gd = A^uj^ [291147]. The inset of Fig. 0 shows 
Ac(S) for the wide container and confirms that this is 
indeed the case. We cannot distinguish any effect of con¬ 
finement, which could be identified as plateaus in the 
increase of Ac; this is to be expected in the Ac <C lx limit, 
which is the case for the wide container. In contrast, in 
the narrow container solid walls fix kc, while with PBC 
the behaviour is not clear, roughly increasing before the 
transition point and then decreasing in a non-monotonic 
way; the uncertainty in the measurements does not al¬ 
low a more accurate conclusion in this case. The marked 
difference between both boundary condition cases sug¬ 
gests that, even though EBC and PBC had equal critical 
points, as measured by (C)^, they do have an influence on 
the modes that are being perturbed. In most of the stud¬ 
ied range k* is consistently higher with PBC, showing 
that solid boundaries can have the originally unexpected 
effect of increasing the critical convection roll size. This 
is due to excluded volume effects near the wall, which de¬ 
crease the density and thus have the effect of exciting a 
lower mode, in our case for A « 25. On the wide container 
wall effects become negligible, and thus k* coincides for 
both types of boundary conditions. 

By taking into account that kc in the wide container is 
not constant, we interpret the LBC transition for lx > A* 
as a series of transitions between energetically similar 
states. Inherent fluctuations are strong enough to allow 
the constant switching between contiguous kc- In terms 
of the relevant scales, this is a conflict between Ac, which 
depends on our control parameter uj, and lx, which is 
fixed. As A* is independent of the container size, the 
critical behaviour for |e| ^ 0 is still expected to be uni¬ 
versal. 

Indeed, Sm(£) shows critical-like behaviour for e < 0, 
as shown in Fig. In the narrow container, both types 
of boundary conditions show the same qualitative growth 
for e < 0, although with consistently lower amplitudes in 
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FIG. 8: Structure factor maximum, Sm, as a function of the 
bifurcation parameter e, for EBC (blue) and PBC (red) in 
narrow (top) and wide (bottom) containers. As reference, 
the best fit for the amplitude of the critical mode is included 
(dashed gray, see main text). 


the EBC case, as previously discussed. For £ > 0 the 
PBC case shows a growth reminiscent of the critical am¬ 
plitude of a supercritical bifurcation. In this case, Sm is 
directly related to the amplitude of the critical mode, as 
the lack of a fixed reference frame makes {n{x,t))^ ho¬ 
mogeneous even in the buoyancy-driven convective state. 
On the contrary, the EBC case immediately decays for 
e > 0. In the wide container both cases coincide within 
error for £ < 0, showing that the discrepancy between 
both cases in the small container is indeed a size-effect. 
Sm again looses significance for £ > 0, and the behaviour 
is erratic due to metastability of the transient region in 
the wide systems. 


E. Dynamics of transient states 

The buoyancy-driven convective state presents com¬ 
plex time evolutions. These are heavily dependent on 
the Ixl^c ratio, based on the constraint that the number 
of convective rolls has to be an integer number. Half¬ 
integer values are possible only with solid wall boundary 
conditions. This implies that non-integer values of Ixl^c 
lead to metastable states, as the number of convective 
cells ric presents intermittent behaviour between the two 
closest values of kc, as also between convection rolls at 
different sides of the container, if walls are present. As 
an example, Fig.[^ shows the temporal evolution of n{x) 
for a system with lx = 80, that is, Ix/^c ~ I-®- ^ 

value of oj just after the transition point, the convective 



FIG. 9: Spatio-temporal contours plot of the number of par¬ 
ticles field, n{x,t), for (a) = 50, (b) = 80 and (c) 

lx ~ 400, with cu = 35 (e ~ 0.06) and EBG. These corre¬ 
spond to Ix/Xc ~ 1, Ix/Xc ~ 1-6 and Ix/X* ~ 8, respectively. 
High density regions are shown in red. Over the middle fig¬ 
ure, the number of convection rolls is indicated for exemplary 
regions. 


cell constantly switches between metastable states; it is 
possible to identify two-rolls and one-roll configurations 
at either side of the system, alternating with no clear pe¬ 
riodicity. We believe this to be an important factor to 
take into account on any study of the dynamics of the 
granular convective state: the size of the container has 
no influence on the critical point of the transition, but 
plays a determining role in the dynamics. In our case, lx 
for the narrow and wide containers was chosen a posteri¬ 
ori to diminish the effects of metastability, considerably 
facilitating the study of precursors. 


As Ixl^c is increased further, a new state becomes pos¬ 
sible at the transition region in which convective cells 
coexist with regions effectively in the Leidenfrost state. 
Fig. shows a period of coexistence, as two pairs of 
convective cells emerge in a confined region of the sys¬ 
tem while the rest remains in the Leidenfrost state. No¬ 
tice how the Leidenfrost region is roughly 200c? wide, far 
larger than A*. We interpret this phenomena as the emer¬ 
gence of a localized state in a non-linear system, a subject 
of increased scientific interest |35j. 






















F. Critical mode amplitude 

It has been shown that both the correlation of the fluc¬ 
tuating velocity field and the critical mode of the density 
fluctuations present critical-like behaviour near the tran¬ 
sition. We now look at the overall transition behaviour 
in the context of bifurcation theory, by following the am¬ 
plitude in the emergent pattern of the critical mode, A^,. 
The emergent pattern is more evident and measured from 
the vertical velocity field Vz{x,t)- Ac is the amplitude of 
the mode kc in Vz{x,t)/^, with kc determined by the 
structure factor maximum. The final value of (^c)t is 
obtained by averaging over the whole simulation time. 

In the seminal work of Swift and Hohenberg, hydro- 
dynamic fluctuations were studied for a molecular fluid 
near the thermal convection instability [49| , and a simple 
model for the Rayleigh-Benard instability was derived. 
In the following we apply the Swift-Hohenberg model to 
the LBC transition, inspired by the evident similarities 
of both phenomena; in terms of bifurcation theory, both 
transitions correspond to spatial-mode selecting bifurca¬ 
tions. Nevertheless, we expect the discrete nature of our 
granular system to have a considerable effect close to the 
transition, manifested as fluctuations arising from the fi¬ 
nite number of particles. Thus, we consider that the 
universal behaviour of the fluctuating vertical velocity 
w{z,t) = Vz{x,t) — {vzix,t))t close to the transition is 
given by the Swift-Hohenberg model for pattern forma¬ 
tion with a stochastic term m. 

dtw = e'w -w^ - {dz;x + klfw + ^/rfC,{x,t), (4) 

with the bifurcation parameter given by e' — kc- In our 
system e' — kc ~ e {as kc ^ 1), and thus in what fol¬ 
lows we take e' = e. Fluctuations are modeled by the 
last term, where ^ is a Gaussian white noise, that is 
{({x,t)C{^'A')) = ~ x')6{t — t'); and r]' is the pa¬ 

rameter of noise intensity m- In our system the zero 
correlation of C is justified by assuming the gaseous phase 
close to the moving plate to be the main source of fluctu¬ 
ations, and to behave strictly as a hard sphere gas. The 
lack of temporal or spatial correlations of the particles 
follows from the the low packing fractions {(j> < 0.2), the 
frequent collisions with the bottom plate compared to the 
mean free flight time, and the randomization of velocities 
due to collision with the dense region. 

It is known that in Q the base state w{x,t) = 0 is 
stable for e' < 0, and presents a supercritical spatial 
instability for e' = 0, which leads to the appearance of 
a pattern, in our case corresponding to convective cells, 
for e' > 0. Following m [50], and confirmed by our 
measured velocity profiles, solutions for the critical mode 
kc can be assumed to be of the form 

w= + + (5) 

V 3 V 3 

with a the amplitude of the pattern with mode kc, depen¬ 
dent on the slow time r = st, and U a general function 




FIG. 10: Amplitude of the critical pattern of the vertical 
velocity field Vz{x,t), Ac, as a function of the bifurcation pa¬ 
rameter e, for EBC (blue) and PBC (red) in narrow (top) and 
wide (bottom) containers. The gray dashed lines correspond 
to fits given by the Swift-Hohenberg model with a stochastic 
term (see main text), with noise level rj = 0.0001. The colored 
dashed lines correspond to fits based on a quintic supercrit¬ 
ical bifurcation, for noise intensity a = 0.0008 in the small 
container systems, and a = 0.001 for the wide cases. 


containing higher order terms in a. Substituting ([^ into 
([^ one reaches the amplitude equation corresponding to 
a stochastic cubic supercritical spatial bifurcation: 

d^a = ea — \a\'^a + y/fj({{T) ( 6 ) 

with r] = Srjk A solution for the probability function of 
a, Ps{\a\,e,r]), can be found from Q and ([^, as shown 
in [13150]. From the shape of Ps the expectation value 
can be obtained |50|, given by 


^max 



(7) 


In our case |aniax| = (^c)f Our measurements are con¬ 
sistent with this form for |e| <C 1, as shown in Fig. 10 for 
narrow and wide systems with PBC and EBC. Neverthe¬ 
less, Q does not capture the shape of {Ac)^ (e) for higher 
values of e, deviating considerably already for e ~ 0.05. 

A higher level of agreement can be obtained by consid¬ 
ering an stochastic quintic supercritical bifurcation (Hj: 


dra = ea— |a|^a -I- ^ r]Vh({{T), (8) 

with h quantifying the strength of the quintic non-linear 
term m- This type of bifurcation may be more rele¬ 
vant for our system, as it has been previously associated 
with parametrically driven spatially extended systems. 
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as Faraday patterns |31j and vertically vibrated series 
of coupled pendula [32]. The former can be considered 
closer to our present system than the Rayleigh-Benard 
scenario, taking into account that the bed of grains is also 
a vertically vibrated medium with a free surface. The 
latter case, on the other hand, could be related to the 
already mentioned low-frequency oscillations [33] . Previ¬ 
ously it was shown that density inverted granular states 
in a quasi-one-dimensional container ly ^ d) be¬ 

have approximately as harmonic oscillators. It can thus 
be inferred that for wider containers, as the ones con¬ 
sidered in this study, the dynamics are analogous to a 
series of coupled oscillators, with a yet unspecified cou¬ 
pling mechanism by shear or other interactions. 

Following a similar method as the previous analysis, 
an expression for the expected value of the amplitude of 
the unstable mode can be obtained (for details of the 
derivation we again refer the reader to [Hi), 

lOmaxI = -f fI/3 (9) 

with fl = (3/4)^/^(9 + ^3(27 — 16/3^))^/^ and f3 = 
with a = rjVh. The shape of © is also shown in 
Fig. now in good agreement for higher e in all cases. 

All systems present the same overall shape of Ac(e), 
with the most significant difference being lower ampli¬ 
tudes for e > 0 in the wide containers. More importantly, 
there is no signihcant difference in the noise intensity 
for all cases, except for the narrow container with EBC, 
where the noise term is lower, cr = 6 x 10“'^ ± 10“®. In 
the narrow container with PBC cr = 9 x 10“^ ± 2 x 10“"*^, 
and cr = 10“^ ± 2 x 10“^ for the wide container with any 
boundary condition. The independence of the noise in¬ 
tensity on N suggests that the relevant quantity for the 
critical dynamics is the amount of particles per critical 
length-scale. A*, which in our cases remains constant. 

IV. CONCLUSIONS 

We have studied the granular Leidenfrost to buoyancy- 
driven convection transition, characterized the precur¬ 
sors and proposed a new interpretation of its universal 
dynamical behaviour. The overall picture is of a contin¬ 
uously fluidized bed of grains which goes from homoge¬ 
neous Leidenfrost conhgurations to increasingly velocity 
correlated convective states, until flows are strong enough 
to sustain the density inhomogeneous buoyancy-driven 
convective state. From a bifurcation theory perspective, 
the convection transition can be understood as a pattern 
formation phase-transition, with the emergence of con¬ 
vective cells with a critical length-scale independent of 
the domain size, which is consistent with previously re¬ 
alized hydrodynamic stability analysis of the Leidenfrost 
state pn] . 

The time-dependent fluctuating convection state can 
be characterized by the correlation time of the fluctu¬ 
ating velocity held, which shows critical-like behaviour 


with an exponent of approximately 0.51. From the self¬ 
correlation it is also possible to observe the inhuence of 
low-frequency oscillations [33] on the huctuating velocity 
held. 

The static structure factor shows the emergence and 
growth of the pattern dominant length-scale. The am¬ 
plitude of the critical mode is also seen to show critical 
behaviour, consistent with a supercritical bifurcation. By 
following the most unstable mode throughout the transi¬ 
tion in wide systems it was possible to conhrm that the 
size of the convective cells is indeed proportional to the 
frequency of energy injection. 

In the transient state of wider systems the Leidenfrost 
and buoyancy-driven convective states can coexist. The 
convective state in this region is constantly evolving, pre¬ 
senting metastability between states with different num¬ 
ber of rolls. As energy increases the stability of the con¬ 
vective cells increases, although their number is deter¬ 
mined by the amount of cells that can be fitted in the 
container. Further increasing the energy leads to a com¬ 
paratively slower process of merging of convective cells. 
The rich dynamics of merging and splitting of convective 
cells in coexistence with the Leidenfrost state in the wide 
systems calls for further research. 

Elastic walls and periodic boundaries present the same 
critical points, disregarding any significant confinement 
effects for containers much bigger than the critical con¬ 
vective wavelength (i.e. larger than 50 particle diameters 
in the cases studied). Slightly dissipative side-walls, on 
the other hand, have the effect of decreasing the amount 
of energy needed to trigger the transition, showing that 
the excitation of the unstable mode at the boundaries has 
a more significant effect than the added dissipation. In 
systems 400 particle diameters wide, in all cases studied, 
the boundary conditions did not have any visible influ¬ 
ence. 

The amplitude of the critical mode of convection is 
seen to be coherent with a quintic supercritical ampli¬ 
tude equation. The agreement is much better than with a 
cubic supercritical bifurcation, associated with the Swift- 
Hohenberg equation. This suggests a new interpretation 
of the transition, closer to spatially extended parametri¬ 
cally driven systems than to Rayleigh-Benard convection. 
We hypothesize that the source of the parametric driv¬ 
ing is not the vibration of the container (which has too 
low amplitude and high frequency to couple with the bed 
dynamics), but the low-frequency oscillations present in 
a density inverted bed of grains, i.e. the granular Lei¬ 
denfrost state. In general, we remark that the universal 
behaviour of the density held can only be captured by 
considering a noise term in the corresponding amplitude 
equation which quantihes the discrete, hnite-number ef¬ 
fects. The noise intensity is seen to be independent on 
the system size, except in the conhned small container. 
This suggests that the transition in wider systems is a lo¬ 
cal phenomenon, with the size of the critical convective 
cell as relevant length-scale. 

As outlook, a derivation of the quintic supercritical 




amplitude equation from a series of coupled oscillators 
with the form derived in |38j would be a way of confirming 
the proposed amplitude equation. 
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